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INTRODUCTION 


In  a  typical  aircraft  gas  turbine  there  are  many  instances  in 
which  rotor  rubs  occur.  TV?o  of  the  most  common  are  blade  tip  and  seal 
rubs,  which  are  caused  by  thermal  mismatch,  rotor  imbalance,  high  "g" 
maneuver  loads,  aerodynamic  forces,  etc.  Current  interest  in  fuel 
efficiency  is  a  consideration  which  drives  the  engine  design  toward 
closer  operating  clearances.  Thus  increasing  the  probablity  of  rotor 
rubs.  The  interaction  of  a  rotor  with  its  case,  (rotor  rubs),  has 
been  studied, to-rof  1  and  3 CL  Ref — \  otudietf  £  steady  state  interaction 
between  a  rotor  with  a  rigid  case  neglecting'friction  at  the  interface 
and  Ba£ — 7  stiidiariC-  a  steady  state  interaction  between  a  linear 
flexible  rotor  and  case  including  friction  at  the  interface^  Ref  1  and  €_ 
2-  did  not  -eoaaidafc  fhe  critical  transient  situation  in  whicli  the  rotor 
bounces  off  the  case.  —  .  I 
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It  is  known  that*  rotor  rubs  can  have  an  important  effect  on  the 
rotor  dynamics.  When  a  rotor  rubs  on  the  case,  a  frictional  force  is 

generated  which  can  drive  a  rotor  to  \rtiirl  in  a  direction  opposite  to 

the  direction  of  rotation,  (backward  whirl).  This  frictional  force  is 
relatively  constant  up  to  the  backward  whirl  speed  at  which  the  rotor 
rolls  around  the  case.  Since  this  rolling  contact  speed  is 
proportional  to  the  rotational  speed  of  the  rotor  times  the  ratio  of 

the  diameter  to  the  rotor  clearance,  the  whirl  speed  can  be  hundreds 

of  times  the  rotational  speed  of  the  rotor;  and  thus  be  potentially 
very  dangerous. 

A\ 

There  are!  two  basic  methods  for  studying  transient  rotor 
dynamics.  One  is  the  modal  method  (ref  3  and  4)  which  expands  the 
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solution  in  terms  of  a  few  of  the  lower  frequency  mode  shapes .  If  the 
transient  under  study  is  localized  (like  a  blade  loss  or  a  rotor  rub), 
the  high  frequency  components  are,  at  least  initially,  dominant.  Thus 
the  modal  method  is  not  applicable  to  this  type  of  transient.  The 
other  method  involves  the  direct  integration  of  the  equations  of 
motion,  which  can  be  done  in  either  of  two  ways,  explicit  or  implicit 
integration.  For  example,  ref  5  used  explicit  integration  of  the 
equation  of  motion,  but  this  solution  is  plagued  with  numerical 
stability  problems.  Further,  ref  6  showed  that  explicit  integration  of 
the  equation  of  motion  was  unstable  when  the  product  of  the  critical 
frequency  (for  any  mode  numerically  possible)  and  the  time  step  was 
large.  Therefore,  the  explicit  integration  can  only  be  done  for  simple 
rotors . 

In  contrast,  the  implicit  integration  tends  to  be  stable  (ref  7 
and  8);  but  it  requires  the  solution  of  a  large  number  of  nonlinear 
simultaneous  equations  at  each  time  step.  Ref  9  used  a  technique 
similar  to  ref  7  except  that  it  was  applied  directly  to  the  second 
order  equation  of  motion.  Ref  9  also  noted  that  the  generalized 
forces  on  a  rotor  were  functions  of  the  generalized  position  and 
velocity  of  the  point  where  the  forces  were  applied  and  its  nearest 
axial  neighbors.  This  allowed  the  variables  to  be  arranged  so  that  the 
Jacobian  of  the  set  of  nonlinear  equations  was  block  tridiagonal. 
Therefore,  computing  time  became  proportional  to  the  number  of 
elements  in  the  rotor  dynamics  model  rather  than  to  the  cube  of  the 
number  of  elements.  The  objective  of  this  study  is  to  refine  the 
method  used  in  ref  9  to  include  an  automatic  time  step  routine;  and 
then  apply  the  technique  to  study  blade  loss  induced  rotor  rubs.  The 
automatic  time  step  routine  is  necessary  so  that  the  time  step  can  be 
varied  as  the  rotor  impacts  the  case.  Also,  the  numerical  stability 
of  the  method  used  in  ref  9  will  be  investigated. 


SYMBOLS 


a  reference  amplitude 

c  radial  clearance 

E  absolute  error  estimate 

F  force 

0  order  of  error  in  Taylor  series 
q  order  of  Taylor  series 

r  radial  displacement 

S  stability  matrix 

t  time 

At  time  step 

u  defined  in  eq.(4) 
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z  independent  variable 

a  given  set  of  constants 

£  damping  ratio 

X  eigenvalue  of  stability  matrix 

p  coefficient  of  friction 

u  frequency 


ANALYSIS 

Numeric al  integration: 

Given  an  arbitrary  vector  function  ^(t)  whose  derivatives  exist 

:pans 
q-k 


4\  V*’ V“  ““  “*■  “  J  ’  ”V  ,'W*  *  —  .v  **«”  ~ 

k  (t)’  a  Taylor  series  expansion  can  be  written: 


2k(t  +  At) 


(At)3 


£  iT~  +  °q-k 


a) 


j-0 


with  remainder  of  order  o  k«  If  the  arbitrary  function  is  chosen  as: 


Z  .  (At)  -Kk) 

Vc  ak! 


(2) 


tho  Taylor  series  for  this  function  becomes: 

q 


K(t  +  At)  ■  E  (i)v°  +  3„ 


(3A) 


j-o 


where  the  binomial  coefficients  are  defined  as 

il 


(0-n 


-  k>'. 


for  j  z  k 
for  j  <  k 


(3B) 


If  the  form  of  the  remainder  is  chosen  as: 


q  "  aku 


(4) 


the  Taylor  aeries  becones; 
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where  the  alphas  are  given  In  ref  7  and  u  can  be  determined  from  the 
equations  of  motion  at  the  advanced  time.  The  form  of  the  set  of  the 
equations  of  motion  at  the  advanced  time  Is : 

£F(r,  r,  r,  t  +  At)  -  0  (6) 


From  the  definition  of  z  ,  the  various  derivatives  become: 


:00 


ak! 


(At) 


k  Zk 


(7) 


Substituting  for  the  various  derivatives  into  the  equations  of  motion; 
and  knowing  the  values  at  the  previous  time,  result  in  the  equations 
of  motion  being  a  function  of: 

£f(u,  t  +  At)  -  0  (8) 

-y 

This  set  of  equations  can  be  solved  for  u  and,  from  this  value  of  u, 
the  remainder  can  be  used  as  an  error  estimate  to  control  the  time 
step.  From  the  definition,  z^  represents  a  nondimensional  form  of  r^. 
Therefore  an  estimate  of  the  maximum  absolute  error  is; 

E-aJIttH  (9) 

where  I  full,  the  vector  norm  is  the  maximum  component  of  u.  The 
computer  code  used  in  ref  9  was  modified  to  include  the  following 
automatic  time  step  althroglm.  If  E>. 01%,  re-do  the  calculation  with 
the  time  step  reduced  by  a  factor  of  10.  If  .01%>E>. 001%,  accept  the 
calculation  but  decrease  the  time  step  by  a  factor  of  2.  If 
.001 %>E>. 0001%,  accept  the  calculation  and  maintain  the  same  time 
step.  If  -0001%>E,  accept  the  calculation  but  increase  the  time  step 
by  a  factor  of  2. 


Numerical  stability: 


The  analysis  of  the  stability  of  the  numerical  integration 
technique  assumes  a  model  of  a  rotor  bearing  system  that  is  linearized 
at  some  instant  of  time.  The  homogeneous  equation  of  motion  for  any 
mode  is: 


if  +  2w£r  +  u>  r  «  0 


(10) 


\rtiere  omega  is  the  natural  frequency  and  zeta  is  the  damping  ratio  for 
the  mode.  For  every  mode  that  is  numerically  possible,  with 
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nonnegative  damping  ratio,  the  amplitude  must  either  remain  constant 
or  decay  in  time*  The  numerical  integration  is  defined  as  unstable  if 
the  amplitude  grows  in  time. 


Substituting  the  Taylor  series  for  the  various  derivatives  into 
the  modal  equation  of  motion  at  the  advanced  time  results  in: 


—I 


fall  -  1)  +  2.1a  At  n+  (m  At)2! 
[_2a2  +  2a jU)  At  5  +  aQ(u>  At)2  J  j 


(t) 


(11) 


j-o 


For  this  value  of  u,  the  Taylor  series  expresses  the  solution  at  the 
advanced  time  in  terms  of  the  solution  at  the  present  time  as; 


**(*♦«>-  Zfe)-- 


j(j  -  1)  +  2ju  At  £+  (ui  At) 


j-0  v 

Defining  the  matrix  S  to  be: 


2a2  +  2a2a>  At  £  +  0 tg(u>  At)' 


Zj(t)  (12) 


/j\  ak[j(j  ~  1)  +  2j(o  At  C+  (m  At) ‘ 
\  /  "  ^2a"2  +  2ajU>At  5  +  a q(u  At)2^J 


(13) 


and  the  vector  Z  whose  kth  element  is  z^,  results  in  the  finite 
difference  equation: 


l(t  +  At)  -  S$(t)  (14) 

This  equation  has  a  solution  of  the  form* 

2(t  +  At)  -  X$(t)  (15) 

where  lambda  is  an  eigenvalue  ofj 

St  -  \t  (16) 

If  the  |X|>1,  the  amplitude  grows  and  the  method  is  numerically 
unstable. 
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Rub  model : 

The  interaction  of  a  rotor  with  its  case  is  a  complicated 
phenomenon.  It  can  involve  non-linear  deformation  of  both  the  rotor 
and  the  case.  Rotor-case  rubs  were  experimentally  studied  in  ref  10. 
Analytically  only  simple  rotor-case  rub  models  are  available; 
therefore*  the  case  was  assumed  to  be  linear  with  dry  friction 
interaction  with  the  rotor.  The  radial  and  tangential  forces  on  the 
rotor  are  then: 


Fr  ■  °»  Fq  -  0  [r(  <  C  (17A) 

Fr  -  -k( | r|  -  C),  FQ  -  yFr  |r|  >  C  (17B) 


RESULTS  AND  DISCUSSION 

The  numerical  method  of  ref  9  employed  a  second  order  integrator 
with  a  constant  time  step.  However,  to  study  blade  loss  induced  rotor 
rubs,  it  is  necessary  to  modify  the  method  of  ref  9  to  include  higher 
order  integrators  with  an  automatic  time  step  routine.  The  automatic 
time  step  routine  is  necessary  so  that  the  time  step  can  be  varied  as 
the  rotor  impacts  the  case.  In  order  to  calculate  high  frequency 
components  accurately,  the  time  step  must  be  less  than  the  period  of 
the  high  frequency  component.  When  only  low  frequency  components  are 
important  the  time  step  can  be  increased  to  decrease  computing  time. 
The  algorithm  discribed  in  the  analysis  section  keeps  the  maximum 
error  in  the  displacement  at  less  than  .01%.  It  tries  to  maintain  the 
error  between  .001%  to  .0001%  by  either  decreasing  or  increasing  the 
time  step. 

Another  way  to  decrease  computing  time  Is  to  use  a  higher  order 
integrator.  Ref  7  studied  the  numerical  stability  of  up  to  a  sixth 
order  integrator  applied  to  a  first  order  differential  equation.  The 
numerical  stability  of  these  integrators  applied  to  a  second  order 
differential  equation  was  given  in  the  analysis  section.  The 
numerical  stability  of  an  integrator  is  based  on  modal  rotor  dynamics 
analysis.  If  the  integrator  is  applied  to  a  mode  which  i~  not  driven 
and  has  damping,  the  amplitude  must  decay  in  time.  Figure  1  shows  a 
stability  map  for  the  integrators  used  in  ref  7  applied  to  a  second 
order  differential  equation.  The  abscissa  is  the  damping  ratio  and  the 
ordinate  is  the  product  of  the  time  step  and  natural  frequency  for  the 
mode.  The  stability  map  has  contours  on  it  for  which  the  amplitude 
does  not  change  from  one  time  to  the  next.  On  one  side  of  the  contour 
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tiie  aiaplitude  grows;  (unstable  region),  and  on  the  other  side  It 
decays,  (stable  region). 

Figure  1  shows  the  stable  regions  for  a  fourth  through  sixth 
order  integrator.  The  second  and  third  order  integrators  were  stable 
everywhere.  For  the  regions  where  the  integrators  were  unstable,  the 
amplitude  grew  by  a  few  percent  per  time  step.  It  would  take  on  the 
order  of  a  hundred  time  steps  for  the  amplitude  to  double,  and  it 
would  take  on  the  order  of  a  thousand  time  steps  for  the  amplitude  to 
increase  by  a  factor  of  a  thousand.  Due  to  round  off  errors,  every 
mode  that  is  numerically  possible  in  the  rotor  dynamics  model,  has  a 
finite  amplitude.  These  amplitudes  may  be  small;  but  if  they  are  in  an 
unstable  region,  in  a  few  thousand  time  steps  they  can  become  very 
large.  For  this  reason,  only  the  second  and  third  order  integrators 
were  used.  This  is  still  a  vast  improvement  over  other  types  of 
integrators  such  as  the  one  used  by  NASTRAN.  NASTRAN  uses  an  implicit 
form  of  the  Newmark-Beta  integrator,  ref  8.  This  integrator  is  second 
order  and  does  not  have  an  error  estimate. 

The  rotor-bearing  system  described  in  ref  11,  (which  dynamically 
simulates  a  typical  small  gas  turbine),  was  used  as  the  example 
problem.  This  rotor  bearing  system  consisted  of  a  shaft  with  three 
disks  mounted  on  two  axially  preloaded  ball  bearings  (fig  2).  In  this 
rotor-bearing  system  the  bearings  were  mounted  in  squeeze-film  damper 
journals,  and  the  journals  had  centering  springs. 

The  first  three  critical  speeds  for  the  rotor  bearing  system 
without  oil  in  the  dampers  are  shown  in  figure  3.  Note  that  all  the 
modes  are  bent-  3haft  modes.  The  "classical"  hierarchy  only  applies 
to  stiff  shafts;  therefore,  the  classical  mode  shapes  do  not 
characterize  the  actual  mode  shapes.  The  first  mode,  about  7600  rpm, 
classically  would  be  the  cylindrical  mode.  But  in  this  case,  it  has  a 
large  amount  of  bending  outward  near  the  shaft  center.  The  second 
mode,  about  9200  rpm,  classically  would  be  the  conical  mode.  In  this 
case,  it  has  a  slight  amount  of  bending  outward  near  the  shaft  ends. 
The  third  mode,  about  11200  rpm,  classically  would  be  the  bending 
mode.  In  this  case,  it  has  a  large  amount  of  bending  throughout  the 
shaft . 

The  rotor-bearing  system  was  modeled  by  using  23  elements.  Prior 
to  the  blade  loss  simulation  the  rotor  was  assumed  to  be  balanced  and 
operating  at  9500  rpm.  The  blade  loss  was  simulated  by  an 
instantaneous  application  of  5  mils  of  mass  excentricity  in  the  far 
disk.  The  equations  of  motion  for  this  system  were  directly 
integrated  by  the  method  used  in  ref  9  with  a  variable  time  step.  The 
output  was  interpolated  to  equal  time  steps;  (100  time  steps  per  shaft 
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rotation),  and  displayed  on  a  CRT,  figure  4.  The  display  showed  an 
oblique  view  of  the  rotor  bearing  system,  with  the  bearing  center  line 
as  the  oblique  axis.  The  transverse  vibration  is  indicated  by  the 
position  of  the  rotor  centerline.  The  scale  of  the  transverse 
vibration  exaggerates  the  amplitude  of  the  vibration.  The  display  on 
the  CRT  was  photographed  at  each  time  step.  These  photographs  were 
then  shown  as  a  motion  picture. 

Figure  5  shows  the  superposition  of  the  first  ten  frames  of  the 
blade  loss  simulation  without  a  rub.  Initially  the  rotor,  the 
bearing,  and  the  mass  center  line  coincided.  After  the  blade  loss,  a 
traveling  wave  starts  at  the  blade  loss  disk  and  travels  down  the 
rotor.  During  the  time  high  frequency  components  are  dominant, 

(because  the  rotor  as  a  whole  is  not  moving).  A  model  analysis  which 
only  uses  the  lower  modes  cannot  discribe  the  motion  during  this  time 
period. 

Figure  6  shows  the  position  of  the  rotor  for  the  first  six 

rotations  of  the  rotor  after  blade  loss  without  a  rub  taking  place. 
During  the  first  rotation,  the  blade  loss  disk  spirals  out.  During 

the  second  rotation,  the  disk  on  the  other  end  of  the  shaft  spirals 

out.  During  the  third  rotation,  the  center  disk  spirals  out.  After 
this  the  envelope  of  the  rotor  positions,  seems  to  oscillate  in  a 
conical  fashion,  with  a  frequency  of  about  1/4  operating  speed.  This 
beating  seems  to  be  at  a  frequency  difference  between  the  operating 
speed  and  the  1st  critical  Speed.  (Ref  12  experimentally  showed  a 
similiar  beat  frequency  between  the  operating  speed  and  the  critical 
speed.)  During  this  time  the  rotor  shape  resembles  the  third  critical, 
except  that  the  bearing  center  line  is  not  in  the  plane  of  the  rotor - 
The  maximum  amplitude  occurred  on  the  blade  loss  disk  on  the  sixth 
rotation  and  on  the  opposite  disk  on  the  fourth  rotation.  The 
conclusion  drawn  from  this  figure  is  that  if  there  is  clearance  space 
down  the  rotor  and  a  rub  occurs,  it  does  not  necessarily  occur  at  the 
blade  loss  disk  first. 

The  rotor-case  rub  was  simulated  by  surrounding  each  disk  with  a 
shroud  that  had  a  2  rail  radial  clearance  and  a  stiffness  of 
100, 000///in.  The  rub  was  induced  by  a  repeat  of  the  blade  loss 
simulation  with  the  clearance  restrain.  Two  rub  simulations  were  run, 
one  with  a  coefficient  of  friction  of  .1  and  the  other  with  a 
coefficient  of  friction  of  .2. 

Figure  7  shows  the  first  6  rotations  of  the  shaft  after  blade 
loss  for  a  coefficient  of  friction  of  .1.  During  the  first  shaft 
rotation  the  blade  loss  disk  spirals  outward  and  bounces  off  the  case 
four  times.  Each  collision  of  the  rotor  with  the  case  sends  out 
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traveling  waves  down  the  rotor.  These  waves  interact  with  each  other 
causing  the  envelope  of  the  rotor  motion  to  be  very  complicated.  On 
the  second  shaft  rotation  both  outboard  disks  bounce  off  the  case  four 
times.  As  the  rotor  continues  to  turn  the  orbit  becomes  more 
circular.  That  is,  the  rotor-case  interaction  becomes  less  of  a 
bouncing  nature  and  more  of  a  continuous  contact.  The  envelope  of 
the  rotor  motion  seems  to  be  oscillating  in  a  conical  nature;  but  both 
outboard  disks  seem  to  remain  in  contact  with  the  case.  The  rotor 
continues  to  whirl  about  the  bearing  centerline  in  the  rotational 
direction  (forward  whirl).  The  frictional  drag  forces  are  not  large 
enough  to  drive  the  rotor  into  backward  whirl. 

Figure  3  shows  the  first  4  rotations  of  the  shaft  after  blade 
loss  for  a  coefficient  of  friction  of  .2.  The  motion  of  the  rotor  on 
the  first  rotation  is  similar  to  the  .1  coefficient  of  friction  case. 
On  the  second  rotation,  the  blade  loss  disk  has  a  very  hard  collision 
with  the  wall,  causing  the  rotor  to  bend  considerably.  On  the  third 
rotation  the  rotor  whirl  direction  changes  from  forward  to  backward 
whirl  and  the  rotor  whirl  begins  to  accelerate  in  the  backward 
direction.  fti  the  fourth  rotation,  the  rotor  motion  becomes  very 
large  and  it  continues  to  grow  on  succeeding  rotations. 

This  exanple  problem  has  shown  that  small  changes  in  the 
coefficient  of  friction,  (from  .1  to  .2)  can  change  a  rotor  response 
to  a  blade  loss  condition  from  a  relatively  safe  response  to  a 
catastrophic  response.  ibr  seal  rii>s  the  coefficient  of  friction  is 
probably  between  .1  to  .  2.  For  blade  tip  rubs,  this  rub  model  is  not 
accurate.  This  type  of  rub  involves  material  removal,  phase  changes, 
and  or  non-elastic  deformations.  If  this  model  were  to  be  used  in  a 
general  manner,  then  the  coefficient  of  friction  would  probably  be 
greater  then  .2. 

In  conclusion,  this  computer  code  allow®  us  to  look  at  blade  loss 
induced  rotor  rd:>s  and  displays  the  rotor  motion  in  a  motion-picture 
format.  A  10-minute,  16-mill  imeter ,  color,  sound  motion-  pic  ture 
supplement  is  available,  on  loan,  that  shows  the  computer  made  motion 
picture  for  the  blade  loss  induced  rotor  rubs. 

SIMMARY  OF  RESULTS  A  VO  CONCLUSIONS 

A  method  for  direct  integration  of  a  rotor  dynamics  svstem 
experiencing  a  blade  loss  indtced  rotor  rub  was  developed.  Tie 
followin''  conclusions  vrere  drawn: 

1.  The  method  was  numerically  stable  for  any  time  step  up  to  a  third 
order  integrator. 
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2.  The  time  step  was  controlled  so  that  the  maximum  error  was  less 
then  .01%  and  the  probable  error  was  between  .001%  to  .0001%. 

3.  For  the  rotor  typical  of  small  gas  turbines  a  small  change  in  the 
coefficient  of  friction,  (from  .1  to  .2),  caused  the  rotor  to  change 
from  forward  to  backward  whirl  and  to  destroy  itself  in  a  few 
rotations . 

This  method  provides  an  analytical  capability  to  study  the 

susceptibility  of  rotors  to  rub  induced  backward  vAiirl  problems. 
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figure  1.  -  Numerical  stability  of  Gear's  integra¬ 
tion  method  applied  to  a  second  order  differ¬ 
ential  equation  for  a  2nd  thru  6th  order  of  in¬ 
tegration.  The  2nd  and  3rd  order  methods  are 
always  stable. 


